/*
 *
 * Medium-scale New Keynesian Model 
 *
 * Cloyne, Dimsdale and Hürtgen (2024) 								
 * "Are Tax Cuts Contractionary at the Zero Lower Bound? 
 *  Evidence from a Century of Data"  		
 *																			
 * James Cloyne, Nicholas Dimsdale and Patrick Hürtgen, 2024
 */


/*
 *
 * Notes:
 *  - This mod-file is based on an example originally provided by Luca Guerrieri
 *      and Matteo Iacoviello provided at https://www.matteoiacoviello.com/research_files/occbin_20140630.zip
 *  - The INEG constraint should theoretically be log_Invest-log(steady_state(Invest))<0, but this will lead
 *      to numerical issues. Instead we allow for a small negative value of <-0.000001
 *
 * Please note that the following copyright notice only applies to this Dynare
 * implementation of the model.
 */

/*
 * Copyright (C) 2021 Dynare Team
 *
 * This is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * It is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * For A copy of the GNU General Public License,
 * see <https://www.gnu.org/licenses/>.
 */

var y           ${y}$ (long_name='output')
    l           ${l}$ (long_name='hours')
    g           ${\hat g}$ (long_name='government spending')
    r           ${r}$ (long_name='interest rate')
    inv         ${inv}$ (long_name='investment')
    w           ${w}$ (long_name='real wage')
    lam         ${lambda}$ (long_name='marginal utility of consumption')
    rk          ${rk}$ (long_name='return on capital')
    q           ${q}$ (long_name='capital price')
    c           ${c}$ (long_name='consumption') 
    mc          ${mc}$ (long_name='real marginal cost') 
    infc        ${infc}$ (long_name='inflation') 
    k           ${k}$ (long_name='capital')     
    kbar        ${kbar}$ (long_name='capital bar')     
    u           ${u}$ (long_name='capital utilitization')     
    a           ${a}$ (long_name='productivity')
    zs          ${zs}$ (long_name='transfers savers')          
    zn          ${zn}$ (long_name='transfers non-savers')                  
    tk          ${tk}$ (long_name='capital tax')          
    tl          ${tl}$ (long_name='labor tax')          
    tv          ${tv}$ (long_name='revenue tax')          
    tyratio     ${tyratio}$ (long_name='tax revenues to output ratio tax')          
    b           ${b}$ (long_name='bonds')          
    rnot        ${rnot}$ (long_name='nom rate')          
    ub          ${ub}$ (long_name='preferences')
    revenue     ${revenue}$ (long_name='revenue')      
    rr          ${rr}$ (long_name='ex-post real rate')
    tt         ${t}$ (long_name='total tax rates')    
;

varexo et ${\varepsilon_t}$ (long_name='tax shock') , eb ${\varepsilonb_t}$ (long_name='pref shock')   
;

parameters 
    beta   
    alpha     
    psi       
    phipi     
    phiy      
    omegap    
    omegaw    
    etap      
    etaw      
    chip      
    chiw      
    delta     
    gamma     
    theta     
    xi               
    gammazs   
    gammazn   
    tl_ss     
    tk_ss     
    tv_ss     
    gammav    
    gammag    
    gammak    
    gammal    
    rhot1       
    tlshare   
    tlshare2  
    rhor      
    adj       
    gshare    
    sb_ss     
    rhoa      
    rhob      
    rhog      
    rhozs     
    rhozn      
    sig_t    
;

beta    =	0.99;
alpha   =   0.33;
delta   =   0.025;
gamma   =   1.5;        
theta   =   0;
xi      =   2;
psi     =   0.8; 
adj     =   6;
omegap  =   0.67;
omegaw  =   0.67;      
etap    =   0.14;
etaw    =   0.14;
chip    =   0.1;
chiw    =   0.1;
phipi   =   1.5;
phiy    =   0.15;
gammag  =   0;
gammak  =   0;
gammal  =   0;
gammazs =   300;
gammazn =   0;
gammav  =   0;
gshare  =   0.108;
sb_ss   =   4*0.963;
tl_ss   =   0.2262; 
tk_ss   =   0.2375;
tv_ss   =   0.089;
rhor    =   0.7;
rhog    =   0.9;   
rhoa    =   0.7;
rhozs   =   0.7;
rhozn   =   0.7; 
rhob    =   0.95;     
rhot1   =   0.9;
tlshare =   0.44; 
tlshare2=   0.31; 
sig_t   =   0.01;
model(linear);

// Implied steady state parameters
# kappap    =   ((1-beta*omegap)*(1-omegap))/(omegap*((1-tv_ss) -omegap*beta*chip*tv_ss +beta*chip));
# r_ss      =   1/beta;
# rk_ss     =   (1-beta*(1-delta))/(beta*(1-tk_ss));
# w_ss      =   ((1/(1+etap))*alpha^alpha*(1-alpha)^(1-alpha)*rk_ss^(-alpha))^(1/(1-alpha));
# kshare    =   (w_ss/rk_ss*alpha/(1-alpha))^(1-alpha);
# lshare    =   (w_ss/rk_ss*alpha/(1-alpha))^(-alpha);
# invshare  =   delta*kshare;
# cshare    =   1 - gshare - invshare;
# L_ss      =   0.25;
# y_ss      =   (1/lshare)*L_ss;
# c_ss      =   cshare*y_ss;  
# zshare    =   1*((1-r_ss)*sb_ss - gshare + tk_ss*rk_ss*kshare + tl_ss*w_ss*lshare + tv_ss);
# znshare   =   zshare;
# zsshare   =   zshare;
# lam_ss    =   (1-theta)^(-gamma)*c_ss^(-gamma)*(1-beta*theta); %internal habits
# alpha_l   = lam_ss*w_ss*(1-tl_ss)/( L_ss^(xi)*(1+etaw));
# kappaw    =   ((1-beta*omegaw)*(1-omegaw))/(omegaw*(1+beta)*(1 + 1*xi*(1+1/etaw)));
# tyratio_ss= tk_ss*rk_ss*kshare + tl_ss*w_ss*lshare + tv_ss;    

lam = -(gamma*(1+beta*theta^2)/((1-theta)*(1-theta*beta)))*c + (gamma*theta/((1-theta)*(1-theta*beta)))*c(-1) + (gamma*theta*beta/((1-theta)*(1-theta*beta)))*c(+1)  + ub/(1-theta*beta) - (theta*beta/(1-theta*beta))*ub(+1); // Internal habits
lam = r + lam(+1) - infc(+1);
rk - 1/(1-tk_ss)*tk = psi/(1-psi)*u;
q  = lam(+1) - lam + beta*(1-tk_ss)*rk_ss*rk(+1) - beta*rk_ss*tk(+1) + beta*(1-delta)*q(+1);
adj*inv(-1) = adj*(1+beta)*inv - q - beta*adj*inv(+1); 
k    = u + kbar(-1);
kbar = (1-delta)*kbar(-1) + delta*inv;
w    = 1/(1+beta)*w(-1) + beta/(1+beta)*w(+1) + chiw/(1+beta)*infc(-1) - (1+beta*chiw)/(1+beta)*infc + beta/(1+beta)*infc(+1) - kappaw*(w - xi*l-  1/(1-tl_ss)*tl + lam - ub);
y    = a + alpha*k + (1-alpha)*l;
rk   = w + l - k;
mc   = alpha*rk + (1-alpha)*w - a;
infc = +(beta-beta*omegap*tv_ss)/((1-tv_ss) - beta*omegap*chip*tv_ss+chip*beta)*infc(+1) + (1-tv_ss)*chip/((1-tv_ss)- beta*omegap*chip*tv_ss +beta*chip)*infc(-1) + kappap*mc+ (1-omegap)/omegap*1/((1-tv_ss) - beta*omegap*chip*tv_ss+chip*beta)*tv- beta*(1-omegap)*1*1/((1-tv_ss) - beta*omegap*chip*tv_ss+chip*beta)*tv(+1);
r_ss*sb_ss*(r(-1) + b(-1) - infc) + g + zsshare*zs + znshare*zn = sb_ss*b + rk_ss*kshare*tk + tk_ss*rk_ss*kshare*(rk + k) + w_ss*lshare*tl + tl_ss*w_ss*lshare*(w+l) + tv + tv_ss*y; 
y        = cshare*c + invshare*inv + g + (1-tk_ss)*rk_ss*kshare*u;
tyratio_ss*revenue = rk_ss*kshare*tk + tk_ss*rk_ss*kshare*(rk + k) + w_ss*lshare*tl + tl_ss*w_ss*lshare*(w+l) + tv+tv_ss*y;
tyratio  = tyratio_ss*(revenue - y);
rnot     = rhor*rnot(-1) + (1-rhor)*(phipi*infc+phiy*y) ;
rr       = r - infc(+1);

g  = rhog*g(-1);
zn = rhozn*zn(-1) - (1-rhozn)*gammazn*(b(-1)-y(-1));
zs = rhozs*zs(-1) - (1-rhozs)*gammazs*(b(-1)-y(-1));
tl = rhot1*tl(-1) + tlshare*sig_t*et;    
tk = rhot1*tk(-1) + (1-tlshare-tlshare2)*sig_t*et;     
tv = rhot1*tv(-1) + (tlshare2)*sig_t*et;  
tt = tl+tk+tv;
a  = rhoa*a(-1);
ub = rhob*ub(-1) + eb; 

[name = 'Nominal Interest Rate (10)', bind='zlb']
r = -(1/beta-1);
[name = 'Nominal Interest Rate (10)', relax='zlb']
r  = rnot;
end;

occbin_constraints;
name 'zlb'; bind rnot<= -(1/beta-1); relax rnot> -(1/beta-1);
end;

steady_state_model;
c=0; 
inv=0;
g=0;
y=0;
l=0; 
k=0;
kbar=0;
u=0;
a=0;
zs=0;
zn=0; 
tl=0;
tk=0;
tv=0;
b=0; 
r=0;
infc=0; 
q=0;
w=0; 
mc=0; 
rk=0;
lam=0; 
ub=0;
tt=0;
rnot=0;
revenue=0;
tyratio=0;
end;

shocks;
    var et = .01^2; 
    var eb = .01^2;
end;

steady;

shocksequence_pref=[-0.2*ones(1,6)];

shocks(surprise);
var eb ; 
periods 1:6 ;
values (shocksequence_pref) ;

var et ; 
periods 6 ;
values -1;

end;

occbin_setup;

shock_mat.matrix_1=options_.occbin.simul.SHOCKS;
occbin_solver(simul_periods=50,simul_check_ahead_periods=100);

oo2_=oo_;
shocks(surprise);

var eb ; 
periods 1:6;
values (shocksequence_pref);

var et ; 
periods 6 ;
values -.0000001 ;
end;

occbin_setup;
shock_mat.matrix_2=options_.occbin.simul.SHOCKS;
occbin_solver(simul_periods=50,simul_check_ahead_periods=100);


